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A  Halo  orbit  about  a  libration  point  of  a  restricted  three-body  system  provides 
additional  opportunities  for  surveillance,  communication,  and  exploratory  missions  in 
lieu  of  the  classical  spacecraft  orbit.  Historically  libration  point  missions  have  focused 
on  Halo  orbits  and  trajectories  about  the  Sun-Earth  System.  This  thesis  will  focus  on 
libration  point  orbit  solutions  in  the  Earth-Moon  system  using  the  restricted  three  body 
equations  of  motion  with  three  low-thrust  control  functions.  These  classical  dynamics 
are  used  to  design  and  optimize  orbital  trajectories  about  stable  and  unstable  libration 
points  of  the  Earth-Moon  system  using  DIDO,  a  dynamic  optimization  software.  The 
solutions  for  the  optimized  performance  are  based  on  a  quadratic  cost  function.  Specific 
constraints  and  bounds  were  placed  on  the  potential  solution  set  in  order  to  ensure  correct 
target  trajectories.  This  approach  revealed  locally  optimal  solutions  for  orbits  about  a 
stable  and  unstable  libration  point. 
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I.  INTRODUCTION 


Libration  points,  also  referred  to  as  Lagrange  points  in  the  literature  [Refs  1-16], 
represent  equilibrium  positions  in  the  restricted  three-body  problem.  Of  the  five  libration 
points,  two  points,  L4  and  L5,  are  stable,  meaning  that  it  is  possible  for  a  spacecraft  to 
remain  stationary  at  that  point  or  orbit  about  it.  The  co-linear  Lagrange  points  LI,  L2, 
and  L3  are  unstable;  yet  provide  a  sensitive  region  of  stability  about  which  a  spacecraft 
may  orbit.  All  points  are  referenced  from  the  bary center  (‘B’)  of  the  system,  which 
defines  the  origin  in  the  reference  frame  and  represents  the  mass  center  of  the  system. 
Figure  1  illustrates  the  positioning  of  the  libration  points  for  the  Sun-Earth  system,  and 


the  convention  that  shall  be  used  to  identify  each  point  throughout  this  thesis.  In  this 
figure,  the  Sun  represents  the  primary  body  of  the  system,  and  the  secondary  body  in  the 
Earth. 

The  most  common  type  of  orbit  about  a  libration  point  is  generally  referred  to  as  a 

Halo  orbit  [Refs  1-12],  and  provides  addition  opportunities  for  surveillance, 

communications  or  exploratory  missions.  Halo  is  not  an  acronym,  the  orbit  is  so  named 
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because  the  orbital  plane  does  not  intersect  the  main  celestial  body  as  a  classical  orbit 
does.  Instead,  the  orbit  resembles  a  Halo  hovering  overhead  as  shown  in  Figure  2.  The 
advantage  of  this  type  of  orbit  over  the  traditional  orbit  is  that  it  generally  provides  a 
continuous  and  uninterrupted  view  of  its  mission  subject. 


•  LI 

Figure  2.  HALO  Orbit  about  L2  Point  of  Earth-Moon  System 

The  purpose  of  this  thesis  is  to  design  an  optimal  Halo  orbit  about  a  libration 
point  of  the  Earth-Moon  system,  using  the  DIDO  optimization  software,  which  is  a 
MATLAB  application  tool.  This  optimal  solution  method  may  be  additionally  applied  to 
any  general  three-body  system,  and  at  the  appropriate  libration  point.  All  libration 
missions  to  date  have  been  in  the  Sun-Earth  System.  This  thesis  will  attempt  to  exploit 
the  Earth  Moon-System  for  a  future  communications  satellite  mission. 

The  design  criteria  or  specifications  for  the  libration  point  orbits  in  this  thesis  are 
based  on  orbital  period,  bounds,  and  constraints  particular  to  the  Earth-Moon  system. 
This  problem  is  scaled  and  non-dimensionalized,  however  different  masses  yield  unique 
mass  ratios  between  the  primary  and  secondary  bodies,  and  alter  the  dynamics  and 
boundary  conditions  of  the  problem  with  respect  to  libration  point  location  and  orbit 
optimization.  Therefore,  the  characteristics  of  the  system  as  well  as  target  orbits  are 
important  in  shaping  the  design  process. 
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II.  BACKGROUND 


In  the  history  of  the  space  program,  there  have  only  been  six  missions  to  libration 
points,  and  all  have  been  in  the  Sun-Earth  system  [Ref  1-2].  The  first  Lagrange  or 
libration  point  mission  was  the  International  Sun-Earth  Explorer-3  (ISSE-3)  [Ref  3] 
launched  in  1978.  ISSE-3  maintained  a  complex  orbit  shown  in  Figure  4,  about  the  LI 
point  to  the  Sun-Earth  system,  where  it  observed  and 
detected  solar  flares  and  cosmic  gamma  ray  bursts.  The 
Halo  orbit  allowed  the  spacecraft  to  make  observations 
over  one  and  a  half  million  kilometers  closer  to  the  Sun 
than  ISEE-1  and  ISEE-2,  which  were  in  Earth  orbits,  and 
demonstrated  the  advantage  and  flexibility  of  Halo  orbit 
missions.  While  the  two  Earth  based  satellites  re¬ 
entered  atmosphere  at  end  of  life,  ISSE-3  was  renamed 
International  Cometary  Explorer  and  sent  to  rendezvous 
with  the  comet,  Giacobini-Zinner  and  flew  through  its 
the  tail  in  1985. 


Figure  3.  ISSE-3  Spacecraft  [From:  Ref  3] 


Figure  4.  ISEE-3  Mission  Trajectory  [From:  Ref  3] 
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Figure  5.  SOHO  Satellite  and  Trajectory  [From:  Ref  4] 


Perhaps  the  most  famous  Halo  orbit  mission  is  the  Solar  and  Heliospheric 
Observatory  (SOHO),  which  was  launched  in  1995  [Ref  4].  Like  its  predecessor  ISSE-3, 
SOHO  also  orbits  the  LI  point  of  the  Sun-Earth  system  and  is  dedicated  to  an  intensive 
and  continuous  study  of  the  star. 


The  most  unique  libration  point  mission  to  date  has  been  WIND,  which  was 
launched  in  1994  as  part  of  the  Global  Geospace  Science  initiative  [Ref  5].  WIND 
investigated  and  studied  plasma,  and  magnetic  field  effects  in  both  ionispheric  and 
magnetospheric  phenomena,  and  made  baseline  observations  in  the  ecliptic  plane  for 


WIND  Extended  Mission 
April  1998 -April  1999 


S17 

8/19/98 


Si 9  Start  2-month 

11/17/98  outer  loop 


S20 

4/1/99 

Enter 

backflip 


XY  Projection  in  Geocentric  Solar  Ecliptic  (GSE)  Coordinates 


Figure  6.  Extended  WIND  Mission  Trajectory  [From:  Ref  5] 
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future  missions.  Its  initial  trajectory  included  multiple  passes  of  the  Moon  before  settling 
into  a  Halo  orbit  about  the  LI  point  in  the  Sun-Earth  System.  Several  months  later  it 
departed  the  L 1  point  for  an  additional  lunar  swing  by  before  it  initiated  a  series  of  petal 
orbits  taking  it  out  of  the  ecliptic. 


Figure  7.  Advanced  Composition  Explorer  (ACE)  [From:  Ref  6] 

In  the  tradition  of  SOHO,  the  Advanced  Composition  Explorer  (ACE),  launched 
in  1997  [Ref  6],  also  orbits  the  Sun-Earth  system  LI  point,  and  obtains  more  specific  and 
detailed  measurements.  The  Microware  Anisotropy  Probe  (MAP)  was  launched  in  2001 
and  marked  the  first  mission  to  the  L2  point  of  the  Sun-Earth  System,  where  it  looks  deep 
in  to  space  to  decipher  the  age,  geometry,  and  size  of  the  universe  without  the 
obstructions  of  the  Earth,  Sun  or  Moon  [Ref  7]. 


Figure  8.  Microwave  Anisotropy  Probe  (MAP)  [From:  Ref  7] 
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The  most  recent  libration  point  mission  is  NASA’s  Genesis,  which  reached  the  LI 
Sun-Earth  point  in  2001  using  a  Lissajous  Orbit  Insertion  (LOI),  which  resembles  a 
figure  eight  trajectory  [Ref  8].  Genesis  is  collecting  actual  specimens  of  solar  wind 
particles  that  it  is  then  returning  to  Earth.  Future  Halo  mission  include  Darwin,  the 
Infrared  Space  Interferometry  Mission  [Ref  9],  which  like  MAP  will  orbit  the  L2  Sun- 
Earth  point  in  search  of  Earth-like  planets  using  six  telescopes.  Darwin  is  not  scheduled 
to  launch  until  2014. 


Figure  9.  Genesis  Lissajous  Trajectory  [From:  Ref  8] 


Figure  10.  Darwin  Telescope  Flotilla  [From:  Ref  9] 
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III.  HALO  ORBIT  PROBLEM  FORMULATION 


A.  COORDINATE  SYSTEM 


1.  Earth-Moon  System 

The  geometry  for  the  restricted  three-body  problem  consists  of  two  coordinate 
systems,  the  synodic  and  the  barycentric  [Ref  13-14].  The  libration  points  in  any  three- 
body  system  exist  in  the  rotating  synodic  (xs,  ys,  z$)  coordinate  system.  The  barycentric 
frame  is  the  inertial  reference  with  respect  to  the  Sun,  and  is  fixed  at  the  barycenter  of  the 
system.  The  subscript  one  identifies  parameters  associated  with  primary  body;  Earth,  and 
the  subscript  two  identifies  parameters  associated  with  the  secondary  body,  which  for  this 
system  is  the  Moon,  shown  below  in  Figure  10. 


ys 


Figure  11.  Earth-Moon  System  Geometry 


2.  Scaling 

The  variable  and  units  in  the  problem  are  naturally  non-dimensionalized.  This 
problem  is  scaled  using  the  variable  //*,  which  should  not  be  confused  with  the 
gravitational  parameter,//  [Ref  13-14].  The  location  of  the  barycenter  for  the  system  is 

m2 

historically  determined  by  the  ratio  //  *  ,  which  is  both  the  mass  ratio  = - ,  and  the 

ml  + m2 
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ratio  used  to  scale  the  distance  between  the  primary  and  second  body  of  the  system  by 
setting  that  distance  =  1.  For  the  Earth-Moon  system  specifically,  ju*=  0.0122,  where  1 
distance  unit  (1  DU)  is  equal  to  the  distance  between  the  Earth  and  the  moon;  384,400 
km. 


y 


<4 


* 


Figure  12.  Mass  and  Distance  Scaling  for  Earth-Moon  System 


X 

■fr¬ 


Mass 

Distance  from  barycenter 

kg 

scaled 

km 

scaled 

Earth 

5.9742  x  1024 

0.9878 

4690 

0.0122 

Moon 

7.3483  x  1022 

0.0122 

379,710 

0.9878 

iable  1.  Mass  and  Distances  for  Earth-Moon  System 


3.  Spacecraft  Reference  and  Control 

The  controls  of  the  spacecraft  are  simply  defined  by  three  thrust  directions  and  are 
referenced  to  the  synodic  system  (Tx,  Ty,  Tz)  as  shown  on  the  following  page  in  Figure 
12.  In  this  figure,  the  vector  R  is  referenced  from  the  origin,  or  the  bary center  of  the 
system.  The  spacecraft  is  also  referenced  from  the  primary  (ri)  and  secondary  fo)  bodies 
of  the  system  for  the  purpose  of  fonnulating  the  spacecraft  dynamics,  which  are  shown  in 
Figure  13  along  with  their  derivation  for  use  in  the  equations  of  motion.  The  thrust  terms 
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represent  the  control  function  of  the  spacecraft  based  on  accelerations  (ax,  ay,  az). 
Accelerations  are  used  in  the  formulating  the  dynamics  in  order  to  simplify  the  problem 
without  the  need  to  select  or  consider  specific  propulsion  ratings  based  on  predicted  mass 
flow  rates  and  ISP  perfonnance. 


B.  EQUATIONS  OF  MOTION 

The  following  equations  are  the  restricted  three  body  equations  of  motion  tailored 
to  the  problem  [Refs  10,13-14],  and  modified  to  include  an  acceleration  term  (ax,  aY,  az) 
to  represent  the  external  force  on  the  system,  which  is  induced  by  the  thrusting  function 
of  the  spacecraft.  The  constant,  ju*  is  the  mass  ratio  of  the  primary  and  secondary 
celestial  bodies  of  the  system  and  is  defined  as  ju*  =  0.0122,  ri  and  r2  are  respectively 
referenced  from  the  primary  and  secondary  bodies  of  the  system  to  the  spacecraft; 


x  -  2y  -  x  =  - 


y  +  2x-y  =  - 


(1  -  ju*)(x  +  ju*)  ju*(x-l  +  jU*) 

3  3 

'i  r2 

(i -M*)y  M*y  .  „ 

..3  .3  +UY 


+  a  x 


eqn(l) 

eqn  (2) 
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eqn (3) 


„  =  _(l-//*)z_//*z+a 


It  is  important  to  specify  the  spacecraft  position  vectors,  ri  and  r2,  with  respect  to 
their  reference  body.  These  vectors  are  different  and  alter  the  dynamics  of  the  problem 
depending  on  whether  the  spacecraft  is  in  the  positive  or  negative  x  quadrant  of  the 
coordinate  system.  This  thesis  focuses  on  solutions  at  the  L2  and  L4  libration  points, 
whose  locations  for  this  problem  are  defined  in  the  positive  x  quadrant.  The  definition  of 
ri  and  r2  is  illustrated  below  in  Figure  13  and  in  equations  (4-5); 


m. 


spacecraft 


i\  =  tJ(x  +  ju)2  +  y2  +z2  eqn  (4) 

r2  =  ^(x-  jLi  +  \)2  +  y2  +z2  eqn  (5) 

C.  LIBRATION  POINTS 

The  actual  equilibrium  points  in  the  system  are  located  in  the  rotating  coordinate 
system  by  setting  the  out  of  plane  velocity  and  acceleration  to  zero  in  the  restricted  three 
body  equations  of  motion  set  [Refs  13-14].  The  thrust  or  acceleration  term  is  also 
dropped  out  in  order  to  find  the  stationary  libration  points  in  the  rotating  frame. 

j  =  (l-^Xx  +  ^)  +  ^*(*-l  +  ^)  eqn  (6) 
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eqn  (7) 


y  = 


(i -t**)y  ,  M*y 


+- 


Q_  (1  ~M*)z  t 


eqn  (8) 


In  order  for  eqn  (8)  to  be  satisfied,  z  must  equal  zero,  therefore  any  equilibrium 
position  in  the  Lagrange  system  must  be  in  the  same  or  orbital  plane  (xy)  as  the  primary 
(mi)  and  secondary  mass  (m2).  Eqn  (7)  can  be  further  simplified  below  by  setting  y=0  in 
eqn  (9). 

0  =  y(l-^-r^  +  ^)  eqn  (9) 

n  r2 


and  then  solving  for  three  co-linear  Lagrange  points  on  the  x  axis  (LI,  L2,  L3),  which  are 
the  three  real  roots  of  eqn  (10). 


^  (1  -ju*)(x  +  ju*)  ju*(x-l  +  ju*)  _ 


eqn (10) 


Substituting  eqns  (4)  and  (5)  into  eqn  (10)  and  simplifying  yields  the  following  equation; 

,  * 


x_az2n+  m- 


(x  +  /u*)  (x-l  +  //*)2 


eqn  (1 1) 


The  solution  to  eqn  (10)  and  the  locations  of  the  three  co-linear  libration  points 
are  obtained  by  first  finding  the  three  real  roots  to  the  Euler  quintic  equations  [Ref  15] 
shown  in  eqn  (12); 


(nij  +m2 )x5  +(3mj  +2m2 )x4  +(3mj  +m2 )x3 -  (m2  +3m3)x2 
-  (2m2+3m3)x  +(m2+m3)=0 


eqn  (12) 


or  as  Vallado  [Ref  14]  expresses  in  three  equations,  eqns  (13-15)  where  mi  is  the  mass  of 
the  primary  body,  m2  is  the  mass  of  the  secondary  body,  and  m3  is  the  mass  of  the 
spacecraft,  which  is  generally  negligible  in  comparison. 

x5  +(3 - /u*)xA  +  (3-2//*)x3  -ju*x2-2ju*x-ju*  =  0  eqn  (13) 

x5  -(3  -  ju*)x4  +(3-2 //*)x3  -ju  *x2  +  2jU  *x-jU*  =  0  eqn  (14) 
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x5  +  (2  -  /u*)x4  +  (1  +  2 ju*)x 3  -  (1  -  ju*)x2  -  2(1  -  ju*)x  -  (1  -  //*)  =  0  eqn  ( 1 5) 

Using  a  numerical  solution  method,  and  substituting  the  mass  values  for  Earth  as 
the  primary  body,  and  the  moon  as  the  secondary,  the  three  real  roots  of  eqn  (12-15)  are 
found  to  be;  (0.8380,  1.1500,  -1.0050)  [Ref  14].  The  specific  nonnalized  x  coordinates 
of  the  libration  points  for  the  Earth-Moon  system  are  then  shown  below; 

Ll=  (  0.8380,  0,  0)  L2=(  1.1500,0,0)  L3=  ( -1.0050,  0,  0) 

Equations  (9)  and  (10)  can  also  be  used  to  find  the  L4  and  L5  Lagrange  points  by 
setting  ri  =  r2  =  1 .  Lagrange  found  the  general  location  of  these  stable  points  based  on 
the  geometry  of  equilateral  triangles  [Ref  13,14,15]  formed  by  the  primary  and  secondary 
bodies  of  the  system  as  shown  in  Figure  14; 


L4  =  (  n~  1/2 , V3 /2 , 0  )  L5  =  ( //  — 1/2 ,  — s/3 /2 , 0 ) 

L4  * 


Figure  15.  L4  and  L5  Libration  Point  Geometry 


For  the  Earth-Moon  System  theses  coordinates  are  defined  in  scaled  units  as 
(0.4879,  0.8660,  0),  and  (0.4879,  -0.8660,  0)  respectively.  Specific  libration  point 
locations  for  the  Earth-Moon  System  in  terms  of  scaled  units  and  actual  kilometers  are 
summarized  in  Table  2  on  the  following  page. 
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Point 

X 

Actual 

Y 

Z 

Scaled 

(km) 

Scaled 

(km) 

Barycenter 

0 

0 

0 

0 

0 

LI 

0.8380 

322,120 

0 

0 

0 

L2 

1.1500 

442,060 

0 

0 

0 

L3 

-1.0050 

386,320 

0 

0 

0 

L4 

0.4879 

187,550 

0.8660 

332,890 

0 

L5 

0.4879 

187,550 

-0.8660 

-332,890 

0 

Table  2.  Earth-Moon  System  Libration  Point  Numerical  Coordinates 
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IV.  OPTIMAL  CONTROL  PROBLEM  AND  ORBIT 

DESIGN 


The  solution  to  any  optimal  control  problem  is  generally  attained  by  “solving  for 
the  state  and  control  histories  of  a  system  subject  to  constraints  while  minimizing  (or 
maximizing)  some  performance  index.”  [Ref  17]  DIDO  is  an  optimization  software 
package  [Ref  18]  that  runs  within  an  existing  MATLAB  program,  it  “employs  a  powerful 
direct  Legendre  pseudospectral  method  that  exploits  the  sparsity  pattern  of  the  discrete 
Jacobian  by  way  of  the  Nonlinear  Programming  solver  SNOPT”  [Ref  19].  After 
formulating  a  general  problem,  a  user  makes  inputs  using  basic  MATLAB  functions  and 
files  according  to  the  appropriate  DIDO  format.  This  fonnat  or  setup  primarily  consists 
of  basic  optimizing  building  blocks  including  dynamics,  constraints,  events,  bounds,  and 
cost  that  make  up  various  sub-files  and  are  mapped  back  to  the  main  solution  file. 

For  simplicity,  a  dual  approach  was  used  to  tackle  this  problem.  First,  an  optimal 
solution  of  an  orbit  about  the  L4  libration  point  was  sought,  since  this  is  a  stable  point 
where  a  solution  is  more  easily  obtained  than  an  unstable  point.  Next,  the  problem  was 
restructured  to  exploit  the  potential  for  trajectories  about  the  unstable  L2  libration  point. 


Figure  16.  Regions  of  xy  Motion  for  the  Earth-Moon  System  [From:  Ref  14] 
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Figure  16  illustrates  the  differences  between  the  zero  motion  regions  about  the  libration 
points  of  the  Earth-Moon  system  in  the  x-y  and  x-z  planes.  In  this  figure,  motion  across 
curves  (C)  of  lesser  value  may  only  be  attained  with  additional  thrust. 


A.  DYNAMICS 

The  restricted  three  body  equations  of  motion  (eqns  1-3)  determine  the  dynamics 
of  the  problem.  These  dynamics  reside  in  an  exclusive  sub-file  that  contains  the  equations 
of  motion.  In  the  dynamic  constraint  t  ,  is  an  independent  variable,  which  is  usually  but 
not  necessarily  time  [Ref  20]. 


rx 

X 

rx 

rY 

y 

rY 

rz 

X 

State  vector  x  = 

I 

rz 

Control  vector  u  = 

ax 

aY 

VX 

vY 

X 

y 

vx 

a7 

vz 

z 

vr 

ax 

X 

vz 

aY 

y 

az  _ 

z 

x(t)  =  f(x(r),u(r),T)  eqn(16) 


B.  EVENT  CONDITIONS 

The  event  conditions  for  the  problem  are  established  by  assigning  values  to  the 
initial  (0)  and  final  (F)  values  of  the  states  or  boundary  conditions.  For  this  problem,  it 
was  not  necessary  to  assign  any  particular  value  to  these  events.  Instead,  it  was  important 
that  the  initial  and  final  events  equal  each  other,  meaning  that  the  final  position  of  the 
spacecraft  match  it’s  starting  position  in  order  to  signify  a  completed  orbit. 

fXO  -  fXF  =  0  Vxo  -  Vxf  =  0 

r yo  -  rY['  =  0  vyo  -  Vyf  =  0 

rzo  -  tzf  =  0  Vzo  -  Vzf  =  0 
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In  order  to  ensure  the  initial  and  final  conditions  are  equal,  the  value  of  each  event 
condition  is  set  to  zero  in  the  main  file  by  setting  the  both  the  upper  and  lower  bounds  of 
the  event  conditions  to  zero. 


C.  GUESSES 

Initial  guesses  are  required  for  the  initial  and  final  conditions  of  the  states, 
controls,  and  time  in  the  DIDO  problem  fonnulation.  The  guess  does  not  necessarily 
need  to  be  feasible,  and  can  be  a  simple  estimate  or  prediction.  However,  in  the  unstable 
libration  point  solutions,  a  reasonable  guess  was  essential  because  of  its  extreme 
sensitivity.  In  this  case,  where  the  user  may  not  be  confident  in  the  reasonability  of  the 
guess,  a  “bootstrapping”  technique  may  be  used  and  is  applied  to  this  problem.  In  this 
process  an  initial  iteration  is  run  using  a  small  number  of  nodes.  This  initial  run  may 
output  a  crude  or  sub-optimal  solution,  but  is  usually  more  reasonable  than  the  guess. 
This  output  is  fed  back  through  the  optimization  process  again,  where  this  initial  solution 
is  used  as  the  guess  for  the  second  iteration.  The  initial  guesses  for  this  problem  are 
scaled  and  defined  in  Table  3  below.  Guesses  for  time  were  based  on  ;rand  2  n ,  which 
are  typical  periods  for  halo  orbits  [Ref  10,12]. 


Guesses 

Stable  Solution 

Unstable  Solution 

Initial  Iteration 

Second  Iteration 

Initial  Iteration 

Second  Iteration 

States 

Initial 

Final 

Initial 

Final 

Initial 

Final 

Initial 

Final 

rx 

0.4879 

0.4879 

0.4883 

0.4883 

1.1500 

1.1500 

1.0505 

1.0505 

rY 

0.8660 

0.8660 

0.8659 

0.8659 

0 

0 

-0.1465 

-0.1465 

rZ 

0 

0 

-0.0018 

-0.0018 

0 

0 

0.0000 

0.0000 

vx 

0 

0 

0.0003 

0.0003 

0 

0 

-0.0191 

-0.0191 

VY 

0 

0 

0.0001 

0.0001 

0 

0 

0.1889 

0.1889 

VZ 

0 

0 

0 

0 

0 

0 

0.0000 

0.0000 

Controls 

ax 

0 

0 

0.0063 

0.0059 

0 

0 

0.0002 

-0.0006 

aY 

0 

0 

0.0063 

0.0060 

0 

0 

0.0002 

0.0022 

az 

0 

0 

0.0058 

0.0058 

0 

0 

0.0006 

0.0006 

Time 

0 

6.2832 

0 

6.2832 

0 

3.1416 

0 

3.6637 

Table  3.  Guesses  for  Optimal  Solution 
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D.  BOUNDARY  CONSTRAINTS 


The  events,  states,  controls  of  the  problem  are  all  assigned  lower  and  upper 
bounds  in  the  main  program  file  in  order  to  specifically  define  the  problem  and  ensure 
feasible  solutions  are  achieved.  As  discussed  in  Section  B,  the  equations  defined  under 
the  event  conditions  were  set  to  equal  zero  such  that  there  was  no  difference  between  the 
initial  and  final  conditions.  All  values  of  the  states,  and  controls,  in  which  the  DIDO 
optimization  software  could  explore  for  a  solution  were  constrained,  so  that  the  scope  of 
the  problem  was  restricted  within  the  vicinity  of  the  desired  solution.  These  constraints 
were  chosen  to  ensure  that  the  output  was  in  fact  an  orbit  about  the  appropriate  libration 
point,  and  did  not  allow  the  spacecraft  to  venture  towards  an  orbit  of  the  Earth,  Moon,  or 
another  libration  point  by  perfonning  an  unnecessary  thrusting  maneuver.  An  example  of 
an  improperly  bounded  problem  is  shown  in  Figure  15  below,  where  the  solution  seeks  a 


Unstable  l_2  Start 


Figure  17.  Unbounded  Trajectory 
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trajectory  about  the  Earth  and  system  barycenter  (B)  after  orbiting  the  moon  although  it 
began  at  the  unstable  L2  libration  point.  Though  not  optimized,  this  trajectory  might 
prove  useful  in  obtaining  a  solution  for  a  low  thrust  transfer  trajectory  from  Earth  to  a 
Halo  orbit  insertion  orbit  about  the  L2  libration  point  of  the  Earth-Moon  system  and  has 
been  previously  presented  in  [Ref  1 1]  and  is  shown  below  in  Figure  18. 


^  Coordinate  (km) 

Figure  18.  Libration  Point  Orbit  Insertion  [From:  Ref  11] 


Constraints  on  the  events,  states,  and  controls  are  expressed  in  eqns  (17-19)  [Ref 
20]  respectively,  and  all  upper  and  lower  level  bounds,  including  time,  are  scaled  and 
listed  in  Tables  4  and  5  on  the  following  page. 
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e,  <e(x(T0),x(zf),T(),Tf)<eu 

eqn ( 1 7) 

x,  <  x(t)  <  xu 

eqn  ( 1 8) 

u,  <  u(t)  <  uu 

eqn  (19) 

Bounds  for  Time  and  Controls 

Time 

Lower 

Upper 

0 

10000 

Control 

Lower  Upper 

ax 

-5.0 

5.0 

aY 

-5.0 

5.0 

az 

-5.0 

5.0 

Table  4.  Time  and  Control  Bounds 


Bounds  for  States  and  Events 


Stable  Solution 

Unstable  Solution 

Stable  and  Unstable  Solution 

States 

Lower 

Upper 

Lower 

Upper 

Events 

Lower 

Upper 

fx 

0.7 

1.1 

1.0 

1.5 

5 

i 

O 

£ 

0 

0 

rY 

0.3 

0.7 

-0.5 

0.5 

rY0~  rYF 

0 

0 

rZ 

-0.5 

0.5 

-0.5 

0.5 

rZ0—  rZF 

0 

0 

Vx 

-10 

10 

-10 

10 

rX0  '  rXF 

0 

0 

VY 

-10 

10 

-10 

10 

rY0—  rYF 

0 

0 

VZ 

-10 

10 

-10 

10 

rZ0  ~  rZF 

0 

0 

Table  5.  State  and  Event  Bounds 


E.  NODES 

The  nodes  represent  markers  or  discrete  points  that  define  the  states  and  controls 
throughout  the  problem.  In  general,  using  a  higher  number  of  nodes  produces  a  more 
accurate  solution  and  takes  longer  computational  time.  Initially,  a  lower  number  of  nodes 
(approximately  100)  was  used  for  the  crude  preliminary  solution  and  was  fed  into  the 
following  iteration  via  the  bootstrap  technique.  For  the  seconds  iteration  a  higher  number 
of  nodes  was  used  (approximately  200)  since  the  guess  was  more  accurate,  and  therefore 
led  to  a  more  smooth  and  precise  solution.  As  the  problem  was  further  explored  and 
refined  a  higher  number  of  nodes  was  used  for  the  initial  and  bootstrapped  solution 
respectively,  which  was  actually  applied  to  for  both  the  stable  and  unstable  solutions. 
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F.  KNOTS 

Knots  are  used  in  DIDO  as  a  part  of  the  optimization  process  and  are  used  where 
there  exists  a  potential  for  discontinuities  in  the  intermediaries  of  the  problem  and 
typically  at  the  end  point  conditions.  In  order  to  satisfy  the  solution  format,  the  location, 
definition,  upper  and  lower  bounds  must  all  be  identified.  The  number  of  nodes  used  in 
obtaining  a  solution  is  also  defined  in  terms  of  these  knots.  For  this  problem,  knot 
locations  were  assigned  to  the  initial  and  final  values  of  time  (to,  tF)  and  were  defined  as 
‘hard.’  Upper  and  lower  knot  bounds  were  also  defined  for  t0  and  tF.  The  value  of  the 
node  knot  number  was  set  to  the  corresponding  number  of  nodes  for  both  the  initial  and 
bootstrap  solution,  as  discussed  in  the  previous  section. 


G.  COST 

The  key  performance  parameter  by  which  the  solution  is  measured  is  prioritized 
by  the  cost  function.  The  minimization  of  a  particular  perfonnance  index  is  given  in  the 
form  of  the  Bolza  cost  function, 


r/ 

•/[*(•)> «(■). fi)U/]  =  E(x(t0), x(r/),r0,r/)+  |  F(x(t),u(t), r),dr  eqn  (20) 

*0 


where  E  is  the  end  point  cost,  and  evaluates  the  cost  function  at  boundary  times  and  F  is 
the  integral  cost  and  is  evaluated  over  the  time  history  of  the  function  [Ref  20]. 
Ultimately  this  function  is  selected  by  the  preference  of  the  user,  but  two  typical  indices 
of  optimality  are  minimum  fuel  and  minimum  time. 

Conserving  or  minimize  fuel  expenditures  is  nominally  a  standard  priority  for  any 
space  mission.  This  is  accomplished  by  minimizing  control  functions  and  thrust 
requirements  within  the  propulsion  budget  of  a  spacecraft.  In  order  to  get  x  independent 
of  the  propulsion  system,  the  following  quadratic  cost  is  used.  It  is  important  to  note  that 
there  are  multiple  solutions  of  optimality  to  a  problem,  and  that  true  optimality  is  only 
defined  by  the  preference  of  the  user  or  customer. 


J  = 


(tf 


T/ 

—  J  (a~  +  a  2  +  a2)d r 


eqn  (2 1 ) 
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H.  PATH  CONSTRAINT 

In  some  optimization  problems  it  is  necessary  to  impose  a  mixed  state  control  in 
seeking  a  solution.  However,  a  path  constraint  was  not  required  for  this  problem  and  was 
therefore  left  unspecified. 
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V.  STABLE  HALO  ORBIT  RESULTS 


A.  STABLE  L4  SOLUTION 

Due  to  the  complexity  of  obtaining  the  solution  for  the  unstable  Lagrange  points 
(LI,  L2,  L3),  an  orbit  solution  set  was  first  found  for  a  stable  Lagrange  point, 
specifically,  L4.  The  initial  solutions  were  propagated  without  a  control  function  in  order 
to  verily  Keplerian  behavior,  and  are  shown  in  following  figures.  Figures  19  and  20 
represent  the  stable  orbit  solutions  at  L4  with  no  control  functions  and  state  boundaries 
imposed.  These  orbits  are  propagated  out  over  a  period  of  100  consecutive  periods  and  it 
demonstrates  how  the  orbit  expands.  Figure  21  shows  a  bounded  solution  with  no  control 
propagated  for  approximately  ten  revolutions. 

1.  Zero  Control  Solutions 


Figure  19.  Unbounded  3D  Solution-  Zero  Control 
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Unbounded  Solution  L4  -  Zero  Control 


-400  - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 

-400  -300  -200  -100  0  100  200  300  400 

scaled  x  coordinate 


Figure  20.  Unbounded  2D  (xy  plane)  Solution-  Zero  Control 


family  of  Unbound**  Solutions  about  L4 


Figure  2 1 .  Family  of  Orbits  for  Bounded  Solution  -  Zero  Control 
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2.  Controlled  Solution 

The  solution  for  the  constrained  and  controlled  orbit  about  the  stable  L4  libration 
point  is  shown  below  in  Figure  22  in  the  xy  plane,  and  again  in  Figure  26  with  respect  to 
the  Moon.  The  orbit  is  also  shown  relative  to  the  position  of  the  moon  in  the  last  figure 
of  this  section.  These  solutions  were  obtained  using  the  quadratic  cost  function,  and  were 
locally  optimal. 


L4  Stable  Orbit 

1 

1  1  1 - 1  1 

linate 

O 

CD 

.  (  \ 

i_ 

O 

o 

o 

1  0.8 

if) 

■  \  ■ 

0.7 

1  1  1  1  1 

0.3  0.4  0.5  0.6  0.7 

scaled  x  coordinate 

Figure  22.  Stable  L4  Libration  Point  Orbit  in  the  xy  Plane 


The  plots  on  the  following  pages  include  the  state  profdes  in  x-y-z,  the  respective 
velocities,  and  control  functions  for  this  particular  solution  all  plotted  against  the 
nonnalized  time,  which  were  the  nodes.  For  this  solution  set,  one  orbit  corresponds  to 
approximately  2  n ,  and  each  plot  reflects  the  periodicity  of  this  time  scale.  It  should  be 
noted  that  the  profiles  in  the  z  coordinate  appear  erratic,  due  to  their  scale,  which  is 
several  orders  of  magnitude  lower  than  the  x,  and  y  coordinates. 
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XYZ  Profiles 


normalized  time 


Figure  23.  XYZ  Profiles 
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Figure  25.  Controls 


Figure  26.  Stable  L4  point  orbit  with  Respect  to  Moon 


27 


B.  VALIDATION 

In  order  to  verify  the  feasibility  of  the  solution,  the  control  solution  generated  by 
the  DIDO  optimizer  is  propagated  through  an  ordinary  differential  equation  solver  using 
the  same  restricted  three  body  equations  of  motion.  The  MATLAB  function  ODE  45, 
with  the  linear  interpolation  of  the  controls  was  used  in  this  case.  Figure  27  shows  the 
comparison  between  the  propagator  solution  shown  in  red  and  the  DIDO  trajectories  in 
blue.  Numerically,  this  difference  in  variation  between  the  solutions  is  on  the  order  of 
zero  to  1 . 1  kilometers. 
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For  this  problem,  it  can  be  shown  that  the  Hamiltonian  plot  is  flat  and  near  zero. 
The  general  Hamiltonian  expression  is  shown  in  eqn  (23)  below  followed  by  the 
Hamiltonian  plot  for  this  specific  solution.  In  eqn  (22)  X  represents  the  Lagrange 
multipliers  or  costates,  which  are  internal  to  the  DIDO  optimization  solution  [Ref  20]. 


H(x,u,X,r )  =  F(x(t),u(t),t)  +  XT  •  f(x(r),u(  r),r)  eqn(22) 


Figure  28.  Hamiltonian  for  Stable  L4  Libration  Point  Orbit  Solution 
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VI.  UNSTABLE  HALO  ORBIT  RESULTS 


A.  UNSTABLE  L2  SOLUTION 

Once  the  stable  solutions  were  attained,  the  orbits  for  the  unstable  points  were 
tackled  with  greater  ease  and  some  success.  Because  of  the  difference  in  location  and 
stability,  the  L2  solution  required  different  boundary  conditions  and  guesses,  but  similar 
constraints.  These  values  were  presented  with  the  orbit  design  in  Section  IV  along  with 
the  stable  solution  values.  As  mentioned  before,  a  reasonable  guess  for  this  solution  was 
necessary  in  order  to  achieve  feasible  results.  Unlike  the  stable  orbit,  it  was  even 
necessary  to  change  the  structure  of  the  guess  to  resemble  an  orbit  in  the  form  of  a  circle 
for  a  feasible  unstable  solution.  Making  a  circular  guess  about  the  unstable  point,  L2, 
encouraged  a  similar  solution  about  the  libration  point.  All  solutions  for  the  unstable 
points  were  found  to  be  locally  optimal  and  had  a  period  that  corresponded  to  n . 
Solutions  in  the  similar  format  presented  in  Section  V  are  shown  below; 


Unstable  L2  Point  Orbit 


scaled  x  coordinate 


Figure  29.  Unstable  L2  Libration  Point  Orbit  in  the  XY  plane 
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XYZ  Profiles 
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Figure  30.  XYZ  Profiles 


Velocity  Profiles 


Figure  3 1 .  XYZ  Velocity  Profile 
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normalized  time 


Figure  32.  Controls  for  Unstable  Orbit  Solution  about  L2  Libration  Point 


Unstable  L2  Point  Orbit 


scaled  x  coordinate 


Figure  33.  Unstable  L2  Orbit  with  Respect  to  Moon 
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B.  VALIDATION 

The  validity  of  the  unstable  point  solution  was  conducted  in  the  same  manner  as 
the  stable  solution.  The  control  solution  generated  by  the  DIDO  optimizer  was 
propagated  through  the  ODE  45  solver,  with  the  linear  interpolation  of  the  controls  was 
used  in  this  case.  Figure  34  shows  the  comparison  between  the  propagator  solution 
shown  in  red  and  the  DIDO  trajectories  shown  in  blue.  The  error  between  the  DIDO 
solution  and  the  propagator  was  comparable  to  the  stable  solution  error.  The  Hamiltonian 
is  shown  on  the  following  page  in  Figure  35.  Again,  this  solution  was  obtained  by  using 
an  initial  circle  guess  solution  of  100  nodes,  followed  by  a  “bootstrapped”  solution  of  200 
nodes. 
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Hamiltonian 


1.5  2 

normalized  time 


Figure  35.  Hamiltonian  for  Unstable  L2  Libration  Point  Solution 
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VII.  CONCLUSIONS  AND  FUTURE  WORK 


Libration  points  provide  additional  locations  for  spacecraft  orbits  with  no 
obstructions  or  interruptions  in  coverage  due  to  eclipse,  which  are  observed  in  traditional 
orbits.  The  design  of  such  orbits  is  particularly  desirable  for  low  thrust  [Ref  21]  vehicles 
since  small  thrust  magnitudes  on  the  order  of  2  to  3  x  10'  m/s“  are  required  to  maintain 
orbit.  Unstable  libration  points  demand  more  stringent  control  functions  than  the  stable 
points,  which  are  more  sensitive  to  slight  deviations  and  perturbations,  would  require 
accelerations  on  the  order  of  10‘4  and  10'5  m/s2. 

Future  work  related  to  this  thesis  might  include  the  incorporation  and 
optimization  of  the  departure  trajectory  from  Earth  orbit  into  a  Halo  orbit  insertion,  which 
would  also  be  required  to  complete  a  mission  to  any  libration  point.  Additionally,  the 
design  of  a  low  thrust  control  system  or  model  might  be  explored  based  on  the 
requirements  dictated  by  the  controls  of  the  DIDO  optimization  software  for  orbits  about 
the  libration  points. 
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